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We study the statistics of thermodynamic quantities in two related systems with quenched disor- 
der: A (l+l)-dimensional planar lattice of elastic lines in a random potential and the 2-dimensional 
random bond dimer model. The first system is examined by a replica-symmetric Bethe ansatz 
(RBA) while the latter is studied numerically by a polynomial algorithm which circumvents slow 
glassy dynamics. We establish a mapping of the two models which allows for a detailed comparison 
of RBA predictions and simulations. Over a wide range of disorder strength, the effective lattice 
stiffness and cumulants of various thermodynamic quantities in both approaches are found to agree 
excellently. Our comparison provides, for the first time, a detailed quantitative confirmation of the 
replica approach and renders the planar line lattice a unique testing ground for concepts in random 
systems. 
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I. INTRODUCTION 



Quenched disorder is common to many condensed mat- 
ter systems. Examples include spin glasses^, elastic 
structures in a random environment 2 or mesoscopic elec- 
tronic systems^. In spite of a large volume of theoretical 
and experimental work on glasses these systems still pose 
interesting challenges. On the theoretical side, much of 
the studies have focused on phase diagrams, the exis- 
tence of phase transitions and their critical behavior. But 
even if the equilibrium phases are known, thermodynamic 
quantities are in general not accessible. Slow glassy equi- 
libration due to the presence of many metastable states 
on large length scales is the main obstacle to the nu- 
merical and experimental study of disordered systems. 
There exist, however, a number of well-tested analytical 
tools to study the equilibrium properties. Common to 
basically all approaches is that they rely on the intro- 
duction of replicas in order to reestablish translational 
invariance. The replicated system is then usually studied 
by a renormalization group (RG) approach or a Gaus- 
sian variational ansatz (GVA)£ since perturbation the- 
ory completely fails. The first method is designed to 
yield the effective pinning potential at large length scales. 
The latter approach aims at constructing a Gaussian trial 
Hamiltonian which describes the glassy phase. The re- 
sults of the RG approaches are, strictly speaking, valid 
only close the upper critical dimension or close to a crit- 
ical point where randomness becomes irrelevant. On the 
other hand, the GVA has to be combined with the con- 
cept of "replica symmetry breaking" which is not gen- 
erally accepted for elastic structures in random medial. 
Although the results of both approaches are similar, the 
different underlying concepts indicate that the present 
general understanding of disordered systems is still in- 
complete. A more general qualitative picture was devel- 
oped in form of the droplet theory for spin glasses^. It 



provides a phenomenological scaling approach to static 
and dynamics properties of the spin-glass ordered phase. 
The properties of this phase are characterized in terms of 
connected clusters (droplets) of coherently flipped spins 
with minimal free energy. It would be desirable to have 
a detailed quantitative test of this theory. Progress on 
numerical approaches as monte carlo simulations at fi- 
nite temperatures were seriously hampered by the slow 
dynamics. Only recently novel power full polynomial al- 
gorithms became available for the study of large systems 
at any temperature^. 

On the experimental side, most of the effort was de- 
voted to spin glasses and pinned vortex systems. Espe- 
cially the latter class comprises a system which is ac- 
cessible to exact both analytical methods and numeri- 
cal algorithms. It is a randomly pinned planar vortex 
lattice that was highlighted by an experimental study 
of magnetic flux lines threading through a thin film of 
the superconductor 2H-NbSe2^- For a certain class of 
2D random systems, including the flux line lattice, new 
promising approaches have been developed. For a planar 
lattice of non-crossing elastic lines pinned by disorder 
a replica Bethe ansatz (RBA) can be employed, yield- 
ing exact results for thermodynamic quantities and their 
cumulantsiSiii*^. For the related^ random bond dimer 
model the partition function can be calculated exactly 
by a polynomial algorithmic^ without the need to run 
slow relaxation dynamics. Therefore, both method are 
perfectly suited to overcome to drawbacks of the above 
mentioned approaches. However, there are also limita- 
tions to the RBA and the dimer simulations. The first 
does not allow to compute correlation functions whereas 
in the latter the choice of parameters of the related line 
lattice model is restricted, e.g., only one particular den- 
sity of lines can be simulated. 

The aim of the present work is to show that the ther- 
modynamics of the two studied model systems provide 
ideal environments for a quantitative test of replica the- 
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ory including replica symmetry breaking and analytical 
continuation. Here the 2D model systems can be consid- 
ered counterparts of exactly solvable ID quantum sys- 
tems which have advanced the understanding of strongly 
correlated systems in general. Our main result is that 
RBA and dimer model simulations agree so well that they 
prove each other to be reliable and thus allow to explore 
many question in detail that had been unaccessible to 
date. 

The rest of the paper is organized as follows. In the 
following section we introduce the line lattice model and 
the dimer model, and explain their mutual relation and 
their connections to other models. In section ITTT1 we re- 
view briefly the replica Bethe ansatz for the line system 
and summarize the results in the known limiting cases. 
We continue in section IIVI with a detailed quantitative 
comparison of the Bethe ansatz predictions for various 
thermodynamic quantities and the corresponding simu- 
lation data for the dimer model. To do so, the Bethe 
ansatz equations are solved numerically outside the va- 
lidity range of the previously studied limiting cases. We 
close with a summary and discussion of our results in 
section 



II. THE MODELS AND THEIR CONNECTIONS 

A. Vortex system 

We consider an ensemble of directed vortex lines con- 
fined to a plane at average distance a = 1/p. The con- 
figurations of a single directed line is characterized by 
its position Xi(z) since overhangs are forbidden. The 
prefercd path of the line results from the competition 
between elastic energy, measured by the line tension g , 
the line interaction in form of a repulsive pair potential 
U (x) that does not allow the lines to cross and pinning 
by a random impurity potential. With a contact repul- 
sion U[x) = c 8{x),c — > oo, ensuring the noncrossing 
condition, the model remains generic*^ and allows for 
the mapping to a discrete dimer model. Quenched disor- 
der couples locally to the vortices via a random potential 
V(r), which we assume to have zero mean and short- 
range correlations 



l/(r)F(r') = A^(r-r'). (1) 

Whenever the disorder correlation length £<j is the small- 
est scale in the problem it can safely be set to zero. The 
total energy can be written as 

H = [ dz E I f ($) 2 + c E - *i) + v ^ z ) 

Throughout the paper, the disorder average will be de- 
noted by ... and the thermal average by ( ... ) . 



B. Dimer model 

The dimer model is defined as follows: Choose a sub- 
set (whose elements are called dimers) of the bonds on 
a square lattice with lattice constant b and linear size 
L such that every of the L 2 — N lattice sites (labelled 
by (ij)) is touched by exactly one of these dimers, see 
Figs -DEI A square lattice rotated by 45 degrees with lat- 
tice constant b/V2 is formed by the centers of the bonds. 
Its 27V sites shall for convenience also be labelled by (ij) 
and it will be clear from the context if the original lat- 
tice or that of the bonds is parametrized. The reduced 
energy of one such complete covering D of N/2 dimers is 
defined by 

where the sum is over all dimers of D. The bond ener- 
gies eij are randomly drawn from a Gaussian distribu- 
tion with zero mean and unit variance ejjejd = fi(ij),(kl)- 
Td is the dimer temperature and measures the strength 
of disorder. The implementation of the polynomial 
algorithmic^ (with exponent ~ 2) on a 32-processor 
cluster allows to compute thermodynamic quantities nu- 
merically exactly - as opposed to, e.g., Monte Carlo sam- 
pling - for sizes up to L = 512 at typically 6000 disorder 
configurations within a CPU time of days. Merely the 
measurement of the specific heat in Section TlV Dl is not 
covered by the polynomial algorithm. Thermal fluctua- 
tions have to be computed from explicitly sampling over 
a representative set out of exp(7VG/7r) possible dimer 
coverings^, again for up to 6000 disorder configurations. 
This gives, however, also reliable results for sizes up to 
L = 256. The typical accuracy of numerical data at the 
given number of disorder samples is ~ 10~ 5 . For further 
details of the algorithm for the simulation of the dimer 
model see, e.g., Ref. 




FIG. 1: Mapping of a dimer (left) to a line configuration 
(right) via XOR-addition of the reference state (middle). 



C. Connections 

These two apparently diverse models are closely re- 
lated. The random bond dimer model can be mapped 
onto an array of lines that interact via a hard-core repul- 
sion, preventing any line crossing. In Fig.[21this mapping 
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and the connections to related models are sketched. From 
the dimer model (A) on top, the discrete lattice version of 
the line model (C) can be reached directly or via the in- 
termediate random solid-on-solid (SOS) model (B). Both 
the SOS model and the discrete lines have their contin- 
uous counterparts - in the bottom line of the sketch - 
that can be treated analytically: The 2-dimensional ran- 
dom sine-Gordon (RSG) model (D) (being equivalent to 
the random-field XY model without vortices) and the 
continuum (l+l)-dimensional elastic lattice of directed 
lines (E). The parameters and observables of the isotropic 
dimer model and its version with a height profile (SOS) 
can be expressed in terms of the anisotropic vortex line 
lattice. In the following we describe the relations be- 
tween the above mentioned models in detail. We start 
the tour through the mapping table of Fig. [21 from the 
dimer model (A) with the dimer temperature Td as the 
only control parameter. The discrete line lattice can be 
reached easily via the map 

(AC) - Add a given dimer pattern and a regular ref- 
erence dimer covering as given in the middle of Figs.m^] 
with an " exclusive or" (XOR) operation. Only if a given 
bond is covered by either the dimer pattern or the ref- 
erence pattern it shall be covered in the resulting line 
configuration, which will be noncrossing lines at average 
density p — 1/(26). A resulting line lattice configuration 
is shown in the right part of Figs. El For the summed 
energies of all covered bonds in the respective configura- 
tions it holds 



isotropic 



anisotropic 



fli({<J) = H d ({e tJ }) 



})■ 



(4) 



on the occupied bonds of the 
= Eij elsewhere. If the origi- 



Here and in the following, the subscript 'd' stands for 
quantities of the dimer model, while line-lattice quanti- 
ties are denoted by the subscript T. {e^ } stands for a 
given distribution of random energies on all of the bonds 
of the dimer model while the set of random energies {e^} 
are defined as = — 
reference pattern and 
nal random bond energies {e^ } are distributed symmet- 
rically with zero mean, so are the bond energies {e'y } 
defining the discrete line lattice model. 

(AB) - A discrete height profile {hij} can be assigned 
to every plaquette of the square lattice in the following 
wajiii: Every bond is given a sign ±1 such that when 
going through rows or columns of plaquettes the sign of 
the crossed bonds alternates. Starting at a given plaque- 
tte with arbitrary height, one moves to the neighboring 
plaquettes and adds to the height +3 times the bond sign 
if the crossed bond is covered by a dimer, or —1 times 
the bond sign if it is not. The resulting numbers define 
a path-independent height profile {hij}, see Fig. for an 
example. 

(BC) - Define a new discrete height profile {Hij} on 
each plaquette by subtracting the height profile {h™ 1 } as- 
sociated with the reference dimer covering, cf. the middle 
part of Fig. 
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A: random bond dimers 


configs.: covering 
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disorder strength: 


1/3? 


B: solid on solid 




configs.: height hi 




disorder strength: 


1/3? 



D: random sine-Gordon 
configs.: height field h(x, z) 
disorder strength: A 
temperature: T 
stiffness: K 



C: discrete line lattice 
configs.: displacement Uij 
disorder strength: l/Tj 
line density: p = \/(2b) 



E: continuous line lattice 
configs.: displacement u(x, z) 
disorder strength: A 
line density: p = 1/(26) 
temperature: T 
elastic moduli: Cn, C44 = T/b 



FIG. 2: Overview of the models introduced in the text, show- 
ing their degrees of freedom and parameters. Top: discrete 
models; bottom: continuous models ; left: isotropic models, 
right: anisotropic models. 
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FIG. 3: Dimer covering (left), reference state (middle) and 
line configuration (right), with height profiles hij (left), h\j 
(middle) and Hij (right). 



The profile {Hij} is then quantized in steps of width 
4. Plaquettes of constant H form domains of average 
width 2b that may be numbered by k. The domains 
are separated by domain walls which can be considered 
as directed lines. Next, we introduce the displacements 
Ufe(i) of these domain walls from their perfectly aligned 
positions where i = 1, . . . , L and j = 2k. As illustrated 
in Fig. the line displacements are determined by the 
height profile, leading to 



Uk{i) 



H 



2k 



Hi t 2k+1 



2k. 



(6) 



Hi 



.ref 



(5) 



It is important to keep in mind that the mapping to a dis- 
placement field makes sense only for a well-defined initial 
configuration of the domain walls implying a fixed den- 
sity. In the grand canonical ensemble induced by the 
sum over dimer configurations, Uk{i) is a good degree of 
freedom only for the configurations corresponding to the 
mean line density p = 1/(26), which, however, carry the 
dominant weight in the thermodynamic limit where the 
density distribution is sharply peaked. 
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FIG. 4: Illustration of displacements Uk(i) in the discrete 
model. 



( CE) - In a continuum formulation each individual 
line can be considered as a directed polymer with energy 
^ J dz {d z Xi(z)} 2 and line tension g. If one includes the 
non-crossing condition imposed by the underlying dimer 
configurations and the coupling to disorder the contin- 
uum model is that introduced in Eq. . In a second step 
this system can be cast into the form of a 2-dimensional 
elastic system. We define a continuous displacement field 
u(r) so that the line positions are Xi{z) = ia + u(ia, z). If 
we include the line interaction in the elastic energy, the 
Hamiltonian for fluctuations on sufficiently large length 
scales can be written as 



— 



A{y(M 2 + y(M 2 +p(rMr)} (7) 



with compression modulus c\\ — aU"(a), tilt modulus 
C44 = g/a and line density p(r) = . 5(x — x.j{z)). The 
correlations of the random potential V(r) are given by 
Eq. ||TJ. The compression modulus is obtained in elastic 
approximation from the line interaction potential U(x), 

d 2 

For the special contact repulsion U(x) = cS(x), c — > 00, 
imposed by the equivalence to the dimer model, the mi- 
croscopic (short scale) en vanishes. However, on larger 
scales (or for more general microscopic interactions) Cn 
assumes a finite value due to entropic contributions. The 
resulting macroscopic en can be obtained from the free 
energy density / of the model of Eq. byi& 



d 2 

C11 = a— [af(a)\ 



(8) 



(BD) - The continuum limit of the random SOS model 
is the random sine-Gordon RSG model with the reduced 
Hamiltonian 8 



ffRSG = J d 2 r{y(V/i) 2 + ^cos[Qh(r)+a 



(r)] 



where K is the reduced surface stiffness, a(r) is a spa- 
tially uncorrelated and uniformly distributed random 
phase, V the coupling strength of disorder and Q — 
2it/A. The proper periodicity Q of the disorder coupling 
can be understood from the periodicity of the related dis- 
crete line lattice, i.e., via the path ACED in the diagram 
of Fig. [3 The latter sequence of mappings also allows to 
conclude about the value of the stiffness K in terms of 
the parameters of the discrete model, see below. 

(ED) - The two continuum models D and E differ by 
their form of elasticity (isotropic vs. anisotropic) and the 
disorder coupling. If we set h(r) — 4u(r)/a and rescale 
the z-axis in model E according to z — > z' — z^c^/cu 
we obtain the isotropic elasticity of model D, cf. Eq. (jj), 
with the stiffness given by 



K = 



a 2 Vc n C44 
16T ' 



(10) 



The disorder energy of Eq. (J7J) can be transformed to 
the form of Eq. J§J) by expanding the local density of the 
lines p(r) in terms of the line displacement u(r) by use of 
Poisson's summation formula, for details see Ref. |2|. The 
relation h(r) — 4u(r)/a between the line displacement 
field and the (dimensionless) height profile is the contin- 
uum version of Eq. JSJl which can be seen with the use of 
Eq. JSJ and the observation that h™* ~ 2j upon coarse 
graining. This relation is in hindsight, knowing that the 
lines have a mean distance of a = 2b, the reason for the 
choice of Q — 2tt/A in model D, Eq. 10. Due to the co- 
ordinate rescaling, the system size A — L X L Z is changed 
when going from model E to D. If we denote the system 
sizes by A e \ and ArsGj respectively, we obtain the ratio 



^RSG 



(11) 



(9) 



Having explained the relations between the models, we 
have to specify how the three continuum model param- 
eters, (i) line stiffness g — CIC44, (ii) disorder strength A 
defined by Eq. and (iii) temperature T, are related 
to the dimer model. First, we observe that only the rela- 
tive strength of the parameters g and A with respect to 
thermal fluctuations is important, i.e., the ratios g/T and 
A/T 2 have to be determined. Out of these we start with 
the reduced single line stiffness g/T. In the dimer model, 
out of the 2L 2 bonds on average likewise L 2 /4 horizontal 
and vertical bonds are occupied. This holds both with 
respect to thermal sampling, i.e., summing over dimer 
configurations, and the disorder average as long as the 
mean random energy is isotropic. For the discrete lines 
after the "XOR" addition with the reference state, an 
average number of L 2 /4 horizontal and L 2 /2 vertical seg- 
ments is implied, corresponding to a mean line density of 
p = 1/(26). We now consider a given discrete line as per- 
forming a one-dimensional random walk with the three 
possibilities of moving to the left or to the right, or to 
rest. The probabilities for the possible steps Xi = — 1, 0, 1 
can be deduced from the average number of occupied hor- 
izontal and vertical bonds. From the above analysis of 
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occupied bonds we find the corresponding probabilities 
Wi = 1/4, 1/2, 1/4 for the steps Xj. For the fluctuations 
of the total horizontal wandering X = J2i x % this leads 
after z/b steps to 



(X 2 



z/b 
i=l 



z 6/2. 



A continuum Hamiltonian for a random walk H 



J dz(d z X) 2 in comparison yields (X 2 ) = -jz and al- 



lows to read off 



9b „ 
— = 2. 
T 



(12) 



The strength of disorder is measured in the continuum 
theory by the variance A of the disorder potential. Not- 
ing that l/Tj is the variance of the (reduced) random 
bond energies in the dimer model we are lead to identify 
\J~K/T with the dimensionless inverse dimer temperature 
l/T^. Allowing for the finite disorder correlation length 
£(2, cf. Eq. which acts as a cutoff in the continuum 
model, to differ slightly from the lattice constant b of the 
dimer model, we have 



A 



kl. 

b Tl 



(13) 



A closer look at the models, however, suggests that this 
relation is valid only if the energies e»j of the horizontal 
bonds of the dimer model are set to zero, i.e., if there is 
no disorder on these bonds. This is because in the con- 
tinuum model of Eq. iJJJ the disorder energy is written as 
the coupling of the local line density p(r) to the disorder 
potential ^(r), 

The random part of the energy is thus not an integral over 
the arc length of the lines but over their z-coordinate. 
Therefore, the disorder energy of a line is proportional to 
its length projected onto the z-axis and not to its over- 
all length. In the original isotropic dimer model, ran- 
dom energy is collected on both vertical and horizontal 
bonds and will thus be proportional to the overall length 
of a fluctuating line. As a consequence, a given disor- 
der in the dimer model corresponds to larger disorder in 
the continuum line model than assumed by Eq. (|13|) . A 
possible method to gauge the disorder strength it to con- 
sider the annealed disorder average —T\nZ of the free 
energy. This quantity can be obtained analytically in 
both models. The calculation for the continuum model 
is straightforward whereas in the dimer model the means 
of Ref. 0| together with the relation of Eq. J2J can be 
used. One finds by comparison between the two models 
the relation 



T 2 b T 2 I 



2GT 2 



2T 2 



e -l/(2T^) 



dx 



arctan x ' 



x 



(14) 



with Catalan's constant G — 0.915966. The contribution 
from the latter two terms in the curly brackets crosses 
over at Td — 1 from zero at small Td to —1/4 at large Td- 
Thus the overall factor of two as compared to the naive 
estimate of Eq. (| 1 31) confirms the conjecture about the 
contribution of the horizontal random bond energies for 
strong disorder (small Td). For weak disorder (large Td) 
the gauging by the annealed free energy yields a factor of 
3/2 compared to Eq. (|13|) . However, in most expressions 
below, the latter two terms will only act as corrections 
at intermediate Td — 1 and are in those cases frequently 
neglected. (At large Td constant additional terms which 
are not related to the gauging procedure dominate.) It 
should be noted that the exact relation Eq. Ijl4|l has been 
derived for one special observable, the annealed free en- 
ergy, and it cannot be expected to hold universally for 
all observables with the same value for the regularization 
length £d- Rather, to different observables the short scale 
modes around the UV cutoff can contribute with a dif- 
ferent weight. The one free parameter £d/6 relating the 
continuum model to the discrete model will therefore be 
considered a fitting parameter, which should, however, 
not turn out to vary dramatically around its expected 
value of order unity from observable to observable. 

Below, results from numerical simulations of both the 
isotropic dimer model with random energies on vertical 
and horizontal bonds and the dimer model with vanish- 
ing energies on the horizontal bonds will be compared to 
theory. For the former, Eq. (|14f) will be used as disor- 
der strength mapping (sometimes without the correction 
terms) while for the latter Eq. I|13fl will prove to fit very 
well. Table[I]summarizes the relations between the dimer 
parameters and the continuum model parameters. 



III. REPLICA BETHE ANSATZ 

The classical statistics of fluctuating elastic lines in d 
dimensions can be described by the quantum statistics of 
a (d — l)-dimensional system of interacting bosons^. In 
the thermodynamic limit, the free energy and its disor- 
der fluctuations are determined by the ground state en- 
ergy of the Bose gas. An important simplification arises 
for (1 + l)-dimensional systems of self-avoiding lines for 
two reasons. First, in the absence of quenched disorder 
the interacting Bose gas is replaced by a 1-dimensional 
free Fermi gas since the non-crossing condition for the 
lines is then automatically fulfilled by the Pauli exclusion 
principl e) 20 ! 21 . Second, 1-dimensional quantum systems 
with sufficiently simple interactions (generated here by 
quenched disorder) can be often treated exactly by Bethe 
ansatz. Indeed, Kardar used the replica method to show 
that the self-avoiding line lattice in a random potential 
maps to a gas of fermions with n spin components with 
SU(n) symmetry interacting via an attractive 5-function 
potential^. The ground state energy of this system can 
be calculated exactly by Bethe ansatz. The analogy be- 
tween the replicated line lattice and SU (n) fermions was 
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model 


line density 


line tension 


disorder 


system size 


continuum line lattice 


p = 1/a 


g/T 


A/r 2 


L x x L z 


isotropic dimer model 


1/(26) 


2/6 




bLx laKT bL 

aq 


anisotropic dimer model 


jj 




- 





TABLE I: Relation between parameters of the continuum line lattice and the dimer model with random vertical and horizontal 
bond energies (isotropic) or only random vertical bond energies (anisotropic). 



examined further in Ref. 1 1 1 L In the following we will 
summarize the main results of the replica Bethe ansatz 
(RBA). 

Upon replication of the system of Eq. J2J) and disorder 
averaging with the aid of Eq. ft}, the equivalent quantum 
system is described by the Eamiltonian 



if- 



_ 2 n N 
" = 2a ^ ^ frr^ 



a</3 j,k 



where N is the number of lines and n the number of 
replicas. In the quantum system the line stiffness g 
corresponds to the fermion mass, the temperature T is 
mapped onto h and the system size in z-direction, L Zl 
onto h/3 where (3 is the inverse quantum temperature. 
We are interested in the ground state with SU{n) sym- 
metry and the corresonding energy Eq (n) . As a function 
of the replica index n, the ground state energy Eq carries 
information on the statistics (cumulants) of the line lat- 
tice free energy. The disorder averaged moments of the 
line lattice partition function are given by 



InZ 



— (-n) j Fi ( 



3=1 



E (n)L z 
T 



n — -ItlZq. (16) 
a 



Here, F J C are the cumulants of the disorder distributed 
line lattice free energy and Zq is the partition function of 
a single line with no disorder contributions. The above 
expansion in the number of replicas n relies on the as- 
sumption of analyticity of the replica free energy (or 
ground state energy) around n — 0. That this assump- 
tion is indeed justified will be demonstrated below. 

The model of Eq. (|15|l is integrable and its ground 
state energy can be obtained in terms of n nested Bethe 
ansatze as demonstrated by Sutherland^ for the repul- 
sive case, and by Takahashi 23 and Kardar— for the at- 
tractive case of interest here. For the Bethe ansatz cal- 
culation n is assumed an integer number. In order to 
extract information from the ground state energy about 
the thermodynamics of the line lattice, Eq(u) has to be 
analytically continued to real valued n. Eq(u) in the limit 
n — > has been obtained in Ref. Using a different 
method for the analytical continuation, an expression for 
Eq (n) at arbitrary n was derived in Ref. ^3- The result 
of the latter reference can be summarized as follows. The 
Bethe ansatz equations can be analytically continued and 
thus yield the allowed wave numbers of the ground state 
wave function for general n. The ground state energy 



can then be expressed in terms of the density function 
g(k) which yields the number L x g{k)dk of allowed wave 
numbers in the interval [k, k + dk]. One finds 



E (n) 



gpA 2 L x ^_ , 



24T 4 



-n(l — n ) + n 



T 2 L. <- K 



'Mi 



dk k 2 g(k). 



K 



The (n dependent) density function g(k) is determined 
by the integral equation 

nk= I dk'g n [{k - k')l d ] g(k') (18a) 

J-K 



with the kernel 



9n[X) 



2 arctan 

m=0 



/?.); 



ni + nm + x 1 



(18b) 



The length = T 3 / (gA) is the characteristic length scale 
of the interaction in the quantum problem. This length 
sets the crossover scale beyond which the line fluctuations 
are dominated by disorder. In the above equations, the 
integral boundary K is fixed by the mean line density via 



K 



p = I dk g(k). 



(18c) 



A" 



In the limit n — > 0, the integral equation l|18af) assumes 
the formiMi 



A" 



1 



K dk> yi d {k-k>) 



- + 7T coth (irl d (k -k'))> g(k') = k 



(19) 

This equation can be solved perturbatively in Kl d , yield- 
ing 



^-Kl d - £r(Kl d ) 3 
2tt 24 v ' 



g(k) = y/1 - (k/K) 2 
144 



{Kl d f 



(20) 



In the opposite limit of vanishing disorder or l d — * oo, 
the exact solution is g(k) — 1/(2tt). Using Eqs. l(T7|l . 
(|18c|) the ground state energy can be calculated from 
Eq. (|20|) perturbatively in Kl d . Interestingly, it can be 
shown order by order in Kl d that the exact result for the 
ground state energy is obtained from the first two terms 



7 



in the square brackets of Eq. (|20() . Indeed one finds thai* 
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K 



K 



dkk 2 g(k) = Y p 3 +p 2 l- 1 (21) 



to any order in pld- Using / = lim n _>o Eo(n)/(nL x ) this 
leads to the mean free energy density of the line system, 



— 7T 2 T 2 3 A 2 

6 g 21 



(22) 



where the contribution from single non-interacting lines 
is given by 



f = C(g)T- 



A 





2£ d T ^ 24T 4 ' 



(23) 



with C(g) a disorder independent contribution that can- 
not be obtained from Bethe ansatz. The result for in- 
teracting lines [Eq. I]22p] was first derived in Ref. 
while the latter result for a single line was first given 
in Refkllo 24] with finite size corrections calculated in 
Refs. |25l26| . It should be noted that in the interaction 
energy of Eq. 1221) the entropic thermal contribution - 
second term - and steric disorder contribution - third 
term - are simply additive with no interference between 
the two effects. They agree in their dependency on the 
system parameters with the respective expressions based 
on scaling arguments 

So far we have only used the integral equation (|18afl in 
the limit n — > 0. However, this equation taken at small 
n contains information about the cumulants of the free 
energy. Therefore, we would like to compute the ground 
state energy E (n) perturbatively innC 1. In general, 
this is not feasible by analytical means. Therefore, below 
we will solve Eq. I|18a|) numerically in order to extract the 
behavior of Eq for small n. There is, however, one impor- 
tant limiting case where analytical progress is possible. 
For strong disorder or low line density with Id -C 1/p, 
the integral Eq. (|18a|) has been shown to assume a par- 
ticular interesting form. It reduces^ after a rescaling 
to the Bethe ansatz equation for the 1-dimensional Bosc 
gas with a repulsive 5-function interactional. The effec- 
tive interaction strength of the Bose gas is proportional 
to n 2 /Id- Therefore, in the interesting limit of small n 
one can either use the numerical solution of the Bethe 
ansatz equations for the Bose gasii or Bogoliubov's per- 
turbation theory 28 to calculate the ground state energy 
Eq (n) of the SU (n) fermions as a polynomial in n to 
lowest order in pld- From Eq(ji) the quenched averaged 
moments of the partition function of the line system can 
be evaluated using Eq. JT5J. For the disorder dependent 
contributions one finds the resulliii 

^=-^(^)\g{^(l-n 2 )ld P + n(ldp) 2 



-^ 2 (M 3/2 +Q-^)^ + 0K)}- 



(24) 



The first term ~ (1 — n 2 )p in the curly brackets describes 
the disorder contribution to the single-line free energy 
cumulants while the following terms stem from line in- 
teractions in the presence of disorder. As can be seen 
from Eq. (|16l) the above result provides the free energy 
cumulants in the dilute limit pld <C 1. At higher densi- 
ties, the mapping to the Bose gas is no longer valid and 
the ground state energy has to calculated by direct nu- 
merical solution of the integral Eq. Q18a[) . This will be 
done in Sec. ITV11 

A natural limit to the validity of the RB A results is set 
by the mapping of the line system with short-ranged cor- 
related disorder to the quantum problem with a singular 
(5-function interaction arising from the disorder average. 
The characteristic length scale in the quantum problem is 
Id which is inversely proportional to the disorder strength 
A. If Id becomes of the order of the cutoff length £<j of the 
short-ranged disorder correlator of Eq. Q , the assump- 
tion of ultra-locally interacting fermions is no longer jus- 
tified. The 5-function interaction, however, is essential 
for the Fermi gas to be solvable by Bethe ansatz. The 
RBA results therefore are valid only for temperatures 



T >T* = (g^A) 



1/3 



(25) 



or, respectively, at sufficiently weak disorder. For lower 
temperatures a modified replica symmetry breaking so- 
lution has been suggested in Ref. [23 for the single line, 
the predictions of which, however, could not be tested to 
date by other means. For the interacting line system the 
modifications at low temperatures implied by the replica 
symmetry breaking solution can be adapted^! However, 
when the RBA results are translated to the random bond 
dimer model it turns out that the low temperature limit 
T < T* can be never realized in the latter model. 



IV. THERMODYNAMICS 

A. Large scale equivalence 

Before studying thermodynamics, we will demonstrate 
that the statistics of the random bond dimer model 
and its associated discrete height profile on large length 
scales are well described by the continuum model for 
elastic lines which can be treated analytically by the 
RBA. From the simulation of the random bond dimer 
model, the correlations of the height profile {hij} of 
model B can be determined very accurately. The quan- 
tity Sh(r) = h(r) — (h(r)) measures the thermal fluctua- 
tions of the height around its state pinned by disorder. A 
Renormalization group calculation predicts for the disor- 
der averaged correlation function 2 ^ 



1 



C(r) = ([Sh(v) - ^(0)] 2 ) = — ln(|r|/6), 



(26) 



i.e., logarithmic growth on large scales, irrespective of 
the value of (h(0)). A possible transition from a glassy 
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low temperature phase with (h(0)) > to a free ther- 
mal phase with (h(0)) = therefore is not reflected di- 
rectly in this correlation function. However, the coeffi- 
cient l/(wK) will be of interest. The stiffness K obtained 
from a measurement of the correlation function for large 
|r| is the large-scale effective stiffness, renormalized by 
contributions from thermal and disorder fluctuations on 
smaller scales. It can be calculated exactly from the RBA 
free energy Eq. (|22|l . First one can use the thermody- 
namic definition of the compression modulus of Eq. (jSJ) , 

en = ad 2 a [a /(a)] , 

and then gets via the relation of Eq. i|l(J|) the effective 
stiffness 



7r / aqA 
A = 16 1 + 



1/2 



(27) 



Note that due to its linear dependence on density, the 
single line free energy of Eq. I|22|) does not contribute 
to K. The result of Eq. ljT7l) is compared in Fig. \5\ to 
the numerical result as obtained from the dimer model 
with isotropic disorder (see also Ref. f8j) and from the 
dimer model with no random energies on the horizontal 
bonds, i.e., e^g/i = 0. The disorder strength A has been 
mapped here to Td according to Eq. (fb^l and Eq. lJT3Jl. 
respectively. The agreement is very good over orders of 
magnitude while the A-T^ relation from the comparison 
of annealed free energy averages seems to fit better than 
the naive estimate. The validity of the replica Bethe 
Ansatz calculation including its sometimes debated ana- 
lytical continuation in n is nicely confirmed by our com- 
parison. Interestingly, no deviation from the RBA result 
of Eq. I|27|l is found in the large disorder, i.e., low Tg 
limit. In terms of the line lattice temperature T modifi- 
cations of the result of Eq. (12211 must occur at T < T* , 
cf. Eq. I|25l) . since otherwise the free energy would not 
converge to a finite ground state energy for T — > 0. How- 
ever, this problem no longer exists after the mapping to 
the dimer model since the ratio V~A/T of disorder and 
thermal energy is controlled by the single parameter T^ 1 
in the dimer model. Thus the pinning strength and ther- 
mal fluctuations cannot be varied independently, and the 
crossover temperature T* can vanish. 

In the pure limit Td — > oo the stiffness approaches the 
value K = n/16 in agreement both with the exact cal- 
culation in terms of the nonrandom dimer modeUi and 
the mapping of the line lattice without disorder to free 
fcrmions. Hence the most accurate simulation of a system 
of noncrossing lines by the dimer model is demonstrated. 
The precise value K = tt/16 is of physical significance as 
shown by the renormalization group (RG) scaling dimen- 
sion of disorder in RSG model Eq. @ 



A A = 2(1 - tt/16 K- 1 ), 



(28) 



see also Ref. The limiting value of K = 7r/16 in- 

dicates that infinitesimal disorder is marginal and the 
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FIG. 5: Large scale stiffness K as a function of the inverse 
disorder strength Td- Shown are the simulation results for the 
dimer model with isotropic random energies and for the model 
with random energies only on the vertical bonds together with 
the RBA result of Eq. Q27p. In the isotropic case Td and A 
are related by Eq. 1141 with £d = 1.15 while in the anisotropic 
case Eq. 1131 with £d = 0.83 has been used. 



system is thus on the borderline between a glassy and 
a thermal free phase. Any finite amount of disorder in- 
creases K which in turn renders disorder a relevant per- 
turbation, leading to a glassy phase. This is consistent 
with the finding of Ref. that the correlation function of 
the height profile, i.e., the one of h and not of 5h, always 
indicates a low temperature glassy behavior. 



B. Free energy, internal energy and entropy 

We now come to a direct comparison of fundamental 
thermodynamic quantities of the dimer model and of the 
continuum model for the line lattice. We start with the 
disorder averaged free energies F = —T\nZ. Due to the 
different meanings of temperature in the line and dimer 
context we will focus on the logarithm of the partition 
functions In Z. When relating the systems we remember 
the energy relation Eq. between dimer and line con- 
figurations. Therefore, the partition functions Zd of the 
dimer model and Z; of the line lattice arc related by 



In Z d + E lci /T d = In Z t 



(29) 



In the disorder average In Zd + E rc f /T d the reference en- 
ergy E le { = H Te f({e' i: j}), cf. Eq. Q, does not contribute 
as the bond energies are drawn from a Gaussian distri- 
bution with zero mean. Higher moments of the free en- 
ergy, however, contain contributions from E ie {. In the 
simulations of the dimer model the moments of the den- 
sity ln(Zj) / (bL) 2 are measured by taking into account 
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the reference energy E le {. The RBA, however, provides 
the statistics of the free energy density of the line lattice 
[Eq. 1)22(1 ] which we denote by fi in the following. Both 
quantities are related by 



In Zi 



Aa T' 



(30) 



where Ai = L X L Z and Ad = (bL) 2 are the system sizes 
of the line lattice and dimer model, respectively. Here 
we have to pay attention to the fact that the two models 
are only equivalent after a rescaling of the z coordinate 
as explained for the mapping between the models D and 
E in Sec. Ill CI According to Eq. l(TT|) we have 



Ai __ a g 
A d ~ 16KT 



a g 

7T T 



1 



agA 

7T 2 T 3 



-1/2 



which has to be used in Eq. i|30|) . The disorder indepen- 
dent part of the single line mean free energy f cannot 
be calculated unambiguously by the RBA. We thus com- 
bine the disorder independent contributions in / ; such 
that in the pure limit (T^ — * oo) the known free energy 
of the dimer model is matched. In this limit, the parti- 
tion functions of the dimer model just counts the number 
of complete dimer coverings of the square lattice. This is 
a complex combinatorial problem as any flip of one dimer 
may necessitate a cascade of flips throughout the system. 
Nevertheless, the result is exactly known to bei& 



lnZd\T d ^oo — —L 2 

IT 



(31) 



in the thermodynamic limit with Catalan's constant G = 
0.915966. Next, we translate the line lattice parameters 
to the dimer model along Table [I] and get for the dimer 
model without horizontal energies, (e^/eh = 0), 



L 2 



1 



G 



1 

T 2 



1 - 



i e d 



a>b 2 Tj 



strength relation between the discrete and the continuum 
model is found, as expected, to be of order one, cf. the 
caption of Fig. 
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FIG. 6: Comparison of the disorder averaged free energy as 
obtained from the dimer model simulations and Eq. I32L 
respectively. In Zi / L 2 is plotted against dimer temperature 
T d ~ T/\TK for the isotropic dimer model (£ d /b = 0.96) 
and the one with random energies only on the vertical bonds 
(£,d/b = 1.00). The simulation data are for system size 
L x L = 256 x 256. 



Why do we have to drop the A 2 -term of the single line 
free energy [Eq. J22Jl] to obtain agreement? One can, in 
fact, imagine a number of reasons for this discrepancy 
between the RBA result and the simulation data. The 
validity of the RBA itself has been critically discussed, 
especially the interchange of thermodynamic limit and 
replica number n — > limit has been questioned. On the 
other hand, the simulations are performed for discrete 
lattice versions of the continuum model which has been 
solved by RBA. We were not able to find a conclusive 



(32) 

If both vertical and horizontal bonds carry random en- 
ergies one should make the approximate replacement 
1/Tj — > 2/Tj, which can be improved by the correc- 
tion terms of Eq. 1|14|) . as explained above. The second 

term ~ T7 2 in the square brackets of Eq. comes answer to what causes the absence of the A 2 -term in the 
from terms proportional to disorder in both the single 
line free energy and the interaction part. The last term 
~ TT comes from the term ~ A 2 in the single line free 
energy. When comparing the result of Eq. (|32|l to the 
simulation results for the dimer model, at first glance 
we find no agreement at all. However, as we will discuss 
shortly, there are indications that one might have to drop 



the TJ~ -term of Eq. l|3^1) . Doing so, we get the plots of 
Fig- ED f° r th e isotropic and the anisotropic random bond 
energies. Only the large T<j limit was fixed by the known 
result of Eq. I|31|l . yet the agreement is excellent over 
orders of magnitude down to small dimer temperatures. 
The only fitting parameter arising from the disorder 



simulation data. But in connection to this it is interest- 
ing to remind of a numerical analysis of the average free 
energy of a single directed polymer in a random potential 
via a transfer matrix method in Ref. l30l In Fig. [7] the 
simulation data of Table II of Ref. |3(3 are plotted, giving 
the average free energy as a function of disorder strength. 
A plot of this kind had not been shown in the cited ref- 
erence. However, it demonstrates that the data obtained 
in Ref. [3C| agree with ours in not finding support for the 
term - A 2 in Eq. (j2l|l . 

Now we compare further thermodynamic quantities for 
the dimer model and the line lattice. The entropy and 
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FIG. 7: Single line average free energy f plotted against dis- 
order strength V 2 ~ A; data taken from Table II of Ref. [30|. 
r ~ 1/g. The solid and dashed curve shows the analytical re- 
sult of Eq. 123H without and with the term ~ A 2 , respectively. 
The reason for the discrepancy is discussed in the text. 



the internal energy of the dimer model 

Sd = 



_ _dT d _ dT c 
U~ d = F d ~ + T d ~S~ d 

are easily be calculated from the RBA. In the dimer 
model simulations, the quenched averaged internal en- 
ergy U d = P{ij) e ij with the disorder configuration 
dependent dimer occupation probability p(ij) of bond 
(ij) can be obtained quite easily since the polynomial al- 
gorithm allows to calculate the probabilities p(ij). The 
entropy is then obtained from the free energy by sub- 
traction. A comparison of the data for internal energy 
and entropy with the RBA prediction is given in Fig. [S] 
where we used the result of Eq. (|32l) and In Z d = In Z; to 
calculate S d and U d . We find excellent agreement with 
£d/b — 1. The slope of S d at T d = is calculated from 
Eq. O to be 



1 d 



L 2 dT d 



S d \ 



G 



1/2 



3/2 



1 



k 

2b, 
(33) 

and matches the simulation data very well. 

Summarizing, the quantitative agreement between the 
RBA results for the line lattice and simulation data for 
the dimer model is very satisfying; it is even more surpris- 
ing for the thermodynamic potentials than for the large 
scale stiffness. The latter is expected to show universality 
in the sense that it does not depend on microscopic de- 
tails of the model, while the former receive contributions 
from all scales. A priori, the sensitivity to the contribu- 
tion from modes close to the UV cutoff might have been 
expected to be important. However, our above results 




FIG. 8: Comparison between simulation data and RBA result 
for the quenched average entropy density Sd/L 2 (£ d /b = 0.98) 
and internal energy density U d /L 2 (£ d /b = 1.00). The data 
are for system size L = 256. 



indicate that the effect of small scales can be simply ac- 
counted for by the single fit parameter t; d /b which was 
found to be very close to the naively expected value of 
one. 



C. Higher moments 

Higher cumulants of thermodynamic quantities de- 
scribe sample-to-sample fluctuations in experimental se- 
tups of mesoscopic dimensions while for macroscopic sys- 
tems their scaling will give information on the selfaverag- 
ing behavior. Analytic expressions for higher cumulants 
of the free energy are available 11 only in the strong dis- 
order or dilute limit l d p <C 1, see Eq. <|24[) . At higher 
densities we have to resort to a numerical computation 
of the ground state energy which then yields In Z n via 
Eq. i|16|) j and thus the cumulants of the free energy. In 
the following we are interested in a polynomial expression 
for E (n). Thus we have to expand in Eq. (|17l) the inte- 
gral of the kinetic energy with respect to n. We introduce 
the dimensionless integral 



E^Sn) = b 6 \ dk k z g(k) = (bl<y / dy y'g(y) 

-K J -1 



^2e j+1 {pl d )n 3 

j>0 



(34) 



defining the expansion coefficients Ej(pld) which depend 
only on the dimensionless parameter pl d . With this def- 
inition, explicit formulas for the variance (second cu- 
mulant) and the skewness (third cumulant) of the re- 
duced free energy can be obtained from the RBA. Using 
Eqs. lfTo|) . ljT7|) together with Eq. and the perturba- 
tive result of Eq. I)24|l we get 



11 



Ai Tb 2 , , , Ail ._ . 

" I (_*m 8/2 ^v /a - — 2 (Wb)3/2 (,,;>l 



A d 3n\aTy \T J A d 3n Tj (perturbative, T rf << 1) 



(mz 3 J T 4a -T £3( ^r ^{l^^ £3(T ^ (exactRBA) 




L2 A l ( 3 l\{A\ 2 bgb A t {3 l\(Ub) 2 , . ... _ n 

11111 11 (perturbative, <C 1) 



A d Vtt 2 4; T 4 



(36) 



Here, the mapping gb/T -» 2, A/T 2 -> (^/fcJ/Tj, 
a — > 2b, (oZd — > (b / £,d)T d / 4 between the line lattice and 
the dimer model with random energies only on the ver- 
tical bonds has been applied. Also the rescaling of the 
volume according to Eq. (|llfl has to be applied for the 
comparison. In terms of the dimer temperature T d , the 
rescaling factor reads 

M = ag f agAy 1 ' 2 4 / 4g d /b\~ 1/2 

A d nT\ 7T 2 T 3 J ^ TT\ 7T 2 T 2 J 

In order to compare the above RBA results to simu- 
lation data for the dimer model over the whole range of 
disorder strength, we solve the Eqs. (|18a|l . (|18c|l for p(k) 
numerically. From p(k) the kinetic energy of Eq. 1|34|) 
and thus the expansion coefficients Ej are obtained. For 
a numerical treatment it is useful to rewrite Eq. (|18a|) in 
the form 

1 f 1 

V = - dy' g n [Kl d (y-y')]g(y') (37) 
n J_x 

with y = k/K and g(y) — g(Ky). At fixed Kl d , the 
dimensionless function g(y) is computed by the inver- 
sion of the discretized integral equation. This inversion 
is quite delicate since the present kind of inverse prob- 
lem - a Fredholm integral equation of the first kind - is 
extremely badly conditioned. An adequate treatment is, 
however, possible by use of, e.g., the method of singu- 
lar value decomposition (SVD) of the discretized integral 
kernel^I. The solution g{y) can only be computed for a 
given value of Kl d . With the so obtained solution we can 
calculate the right hand side of 

J7 = J dy g(y). 

The dimer model density p — 1/(26) then fixes K and 
we thus obtain the disorder strength l/l d that had been 
implied by our initially chosen value for Kl d . In this 
approach we cannot - due to the coupling of the BA 
equations - modify n and l d independently. A modi- 
fied n implies a modified K, which results in a different 



value of l d . We hence have to adjust the parameter Kl d 
in Eq. i|37f) upon change of n such that l d remains con- 
stant. In practice, this is realized by the simple method of 
nested intervals. The necessary number of discretization 
points for the integral Eq. I|37l) depends crucially on con- 
sidered range of parameters. An increase in the length 
scale l d stretches the integral kernel g n (kl d )/n along the 
abscissa, a decreasing n does so along the ordinate. Thus, 
in order to keep up a given level of accuracy, the number 
of discretization points has to increase like l d /n. We are 
unfortunately interested in small n as we want to extract 
the coefficients Sj from the behavior around n = and 
moreover in large l d as the result for small l d is known 
analytically. Therefore it is important that the reliability 
of the numerics can be checked in limiting cases where 
analytical results for g(k) are available. The limit of large 
l d — * 00 or n — > 1 corresponds to lines without quenched 
disorder which are described by free fermions with a con- 
stant density g(k) — 1/(2tv). In the inset of Fig. [51 a plot 
of the numerical solution in the latter case is displayed, 
showing very good agreement with the analytical expec- 
tation. An other limit which can be compared to ana- 
lytical results corresponds to n — > 0. Then g{k) can be 
calculated perturbatively 11 in Kl d , see Eq. (|20ll . We com- 
pare the numerical solution for n = 10 -3 and Kl d = 0.1 
with the result of Eq. I|2(JI) in Fig. |5J Again the compar- 
ison is satisfactory while the mismatch at small k/K is 
due to the smallness of n which necessitates a high dis- 
cretization level. In practice, n is chosen as not to require 
more than 10 3 discretization points for a relative accu- 
racy of 10 -4 in Eki n (n). Checks against the exact results 
for the quenched averaged free energy [Eq. I|22|l ] and the 
strong disorder limit of the cumulants [Eq. I|24|l ] are also 
satisfactory. 

Having determined the solution p(k), the coefficients 
Ej of Eq. (|34|) can be extracted from the numerically cal- 
culated Eki u (n) by repeated extrapolation to n = 0, sub- 
sequent subtraction of this value from the nnite-n result 
and final division by n. This straightforward procedure 
is the best we could think of but it is still error-prone. 
While the desire for a small extrapolation error requires 
having data points as close as possible to n = 0, a simple 
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FIG. 9: Numerical solution g(k/K) for small ldK, n — 
and for l^K — > oo (inset). The mismatch at small k/K be- 
tween the analytical expression and the numerical solution for 
ldK = 0.1 is reduced with an increasing number of discretiza- 
tion points. For the extraction of cumulants, larger values of 
n/ld have been used where a discretization level of 10 3 points 
suffices. 



calculation of error propagation shows that the error of 
a given data point at finite n scales as n~~ where 
k stands for the order of the coefficient and <5i?km for 
the original error of the data. However, we were able to 
achieve sufficient high accuracy to extract reliable val- 
ues for the second and third cumulant. The above de- 
scribed scheme was implemented with the use of Numer- 
ical Recipe routines^. The data for the smallest dis- 
order values needed the highest discretization level and 
consumed about 200h computation time on a 2GHz pro- 
cessor for the determination of the cumulants at a given 
value for Id. 

Before we compare the RBA predictions for the cu- 
mulants to the simulation results for the dimer model 
we would like to make some remarks on the computa- 
tion of cumulants in the dimer model simulations. In the 
replica theory the cumulants of random thermodynamic 
quantities appear naturally, while in the simulations the 
moments are immediately accessible. Both are related 
as follows. If the generating functional of the moments 
{m p } is 

2 m 2 3 m 3 

M(n) = 1 + nmi + n + n + • ■ ■ , 



then 



lnM(n) = nni + n 2 ^- + n 3 ^- + 



generates the cumulants {f%>}- The lowest cumulants ex- 
pressed in terms of moments are 



«3 



Ml 

m 2 



m3 — 3mi?7i2 + 2m\. 



(38) 



From the dependence of Eq. (|16fl on the system size Ai = 
L x x L z it follows immediately that all the free energy 
cumulants scale linearly in Ai. Apart from the average 
free energy, the reduced cumulants k p /k\ = F% /F P hence 
vanish in the thermodynamic limit as F£ / F P ~ A^~ p as 
one would expect from the central limit theorem. The 
distribution of the free energy becomes infinitely sharp 
in the limit of large systems. In other words, the vortex- 
line array is self averaging which in light of the infinite 
correlation length reflected by logarithmic correlations 
had not been evident a priori. 

For the determination of the cumulants from the sim- 
ulation data the following problem is entailed. The mo- 
ments m p scale with the system size like m p ~ Af . From 
Eq. (|3*5|) we see that a cumulant of order p has to be calcu- 
lated as a sum of terms that grow by a factor A p ~ x faster 
with the system size then the cumulant itself. Therefore, 
at a given accuracy of the simulation data for the m p , 
which primarily depends upon the number of disorder 
samples, a limit is set to the system size up to where 
cumulants can reliably be obtained. This maximum sys- 
tem size decreases with the order of the cumulant. On 
the other hand, finite size effects have to be minimized as 
well and, as a consequence, at the achieved precicion of 
10~ 5 for the moments m p in the dimer model, the vari- 
ance can be trusted only for system up to size L = 64 
and the third cumulant up to size L = 16. 

In the following, we compare the simulation data for 
the cumulants (lnZ ; )£/L 2 for p = 2, 3 with the RBA 
predictions of Eqs. (|35|l . I|36|l . The variance (p = 2) as 
a function of Td has been computed for the dimer model 
of size L = 64 and is shown in Figs. EH IHI together 
with the RBA result. With the only fitting parameter 
£d/b = 0.8 we find very good agreement. In addition, the 
plots show that our numerical solution of the full Bethe 
ansatz equations nicely confirms the perturbative solu- 
tion at small Td. The deviation at larger Td shows that 
the numerical solution is inevitable for a comparison with 
the simulation data. A closer look at the analytical re- 
sult of Eq. (|35|) and the data points of Figs. 1101 II ll shows 
that (In Zi) 2 ./ L 2 does not follow strictly a power law in 
Td- The third cumulant (skewness, p = 3) is shown in 
Fig. El For the reasons explained above, the dimer data 
are to be trusted only for the small system L = 16 and 
here only for T<j < 1. In this range, no substantial devi- 
ations from the perturbative evaluation of the cumulant, 
cf. Eq. I|36(l . is expected. Indeed, the agreement between 
theory and simulation is again very good with £d/b = 0.8. 
Although the numerically determined coefficient £3 has 
an error of only ~ 15%, the actual uncertainty of the 
third cumulant is larger. The reason is the following. 
In the exact RBA result of Eq. (|36|) the magnitude of 
the second term ~ £3 amounts throughout the studied 
parameter range to about 85% of the first term. Since 
the terms are subtracted, the original error of 15% gets 
amplified to ~ 100% in the final expression for the third 
cumulant. However, this does not restrict the comparison 
since the perturbative RBA solution could not have been 
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FIG. 10: Simulation data for the dimer model with random 
energies on the vertical bonds only (diamonds), the numerical 
result of the RBA equations (full line) and the perturbative 
RBA result (dashed line) for (In Z t )Z) / L 2 . The fitting param- 
eter is = 0.8, dimer system size is L — 64. 
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FIG. 12: Simulation data (diamonds) and perturbative RBA 
result (full line) for the third cumulant (skewness) (In Zi)%/L 2 . 
The dimer system size is L — 16. Note that the data do not 
follow a power law. 
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FIG. 11: Same as Fig. ^|with a smaller range of T d . 



corrected noticeably in the range T d < 1 which is set by 
the simulation data at hand. Note that both the simula- 
tion data and the perturbative RBA result consistently 
do not obey a power law in Td, see Fig. 1121 

It must be noted that there is agreement for the third 
cumulant only if the term ~ n 3 from the single line con- 
tribution ~ ^(1 — n 2 )A 2 in Eq. is taken serious, 
otherwise not even the sign of the cumulant would match. 
However, the part linear in n stemming from the same 
contribution we had to drop in the comparison of the 
average free energy, being consistent with another in- 
dependent study^ii of the single line free energy. Since 



a single line (or directed polymer) is expected to have 
an asymmetric disorder free energy distribution, we do 
not argue for a complete irrelevance of the contributions 
~ A 2 in the single line terms of Eq. (|24(l . consistent with 
our simulation data for the third cumulant, but suggest 
a reassessment of the free energy distribution of a single 
line in random media. 



D. Specific heat 

The specific heat of disordered systems is strongly in- 
fluenced by the complex nature of thermal excitations 
about the pinned ground state. In the context of spin 
glasses the excitations are considered as droplets, i.e., 
connected regions in which the thermal activated config- 
urations differ from the ground state configuration^*^. 
Droplets appear on all length scales with the lowest en- 
ergy ones appearing on largest length scales. Here we 
would like to test the RBA prediction for the mean spe- 
cific heat by comparing to our simulation results. The 
disorder averaged specific heat of the dimer model with 
random energy on horizontal and vertical bonds can be 
measured via the thermal fluctuations of the dimer en- 
ergy, 



_ T -2 (Hj)-(H d ) 2 

Cd = L, =2 > ( 39 i 

where H d is the Hamiltonian of the dimer model, see 
Sec. Ill Bl From the RBA result for the free energy 
Eq. 12 11 (again without the A 2 -term in the single line 
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free energy) the specific heat is calculated easily using 

^=§^(r d h7^), (40) 
i 



In Fig. ^5] we compare this RBA result to the simulation 
data which we obtained via Eq. I|39|) . We find rather nice 
agreement over the entire range of Td with the choice 
£d/b = 0.98 for the only fitting parameter, which is con- 
sistent with our findings above. There are no further ad- 
justable parameters in the comparison shown in Fig. 1131 
Also note that the dimer specific heat probes for the 
agreement of the RBA and simulations predominantly in 
the region around Td—1, the drop to zero for small and 
large Td being generic rather than specific. So it can be 
considered complementary to the mean free energy which 
tested for amplitude and exponent at small Td while the 
large Td (or pure) limit was fixed to the exactly known 
result, see Section llVBI 
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FIG. 13: Mean dimer specific heat cj- Shown are the simu- 
lation data (diamonds) and the analytical result of Eq. I|41|l 
derived from the RBA (full line). The simulated dimer model 
has random energies on horizontal and vertical bonds, and is 
of size L = 256. 

The linear low temperature behaviour of the specific is 
typical of random systems. It can easily be understood by 
considering the lowest lying excitations on a given length 
scale t. On each scale for a given disorder environment 
the ground state (with zero energy) and the lowest exci- 
tation (with energy E) form a two level system, whose 



where In Zd = In Zi is given by Eq. (|30|l . After the map- 
ping A/T 2 2(£ d /6)/T 2 , bg/T -> 2, p -> 1/(26), the 
mean specific heat of the dimer model reads 



(41) 

I 

specific heat is 

_ d Ee- E I T 
° LE ~ dfe- E l T + l 

The excitation energies obey a length scale dependent 
disorder distribution pi(E) and the contribution to the 
mean specific heat from each scale is 

poo 

5ct= / dE Pl {E)ci :E ~ Pl {Q+)T. 
Jo 

For a finite density of excitations at small energies, 
liniE^o+ pe{E) > 0, the specific heat as a superposition 
of exponentials contributions from each scale will be lin- 
ear at small temperatures which is a famous insight of 
Anderson et alM> Integration over all length scales al- 
lows to write the mean specific heat as c = dt See 
which becomes exact in the limit T — > 0. In the droplet 
theory of spin glasses the distribution function pi{E) is a 
central quantity. From the finite size scaling of the mean 
specific heat one can hope to obtain information on this 
distribution function as with growing system size larger 
droplets will fit into the system. However, the droplets in 
our dimer model simulation seem to be dominated by the 
system boundaries. The number of configurations of low- 
est excitation energy that differ only on the boundaries 
from the ground state configuration scales like the linear 
system size L, yielding a 1/L- decay for the specific heat 
(instead of a growth from bulk droplets) to its asymp- 
totic value . In fact, this scaling behavior is observed for 
the specific heat data of our dimer model simulation, see 
FigO Indeed, low-lying excitations on the boundary 
can be easily identified. Consider a bond on the bound- 
ary that is occupied by a dimer in the ground state. A 
configuration that does not cover this very bond may re- 
main unchanged on all the other bonds since simulations 
are done with open boundary conditions. The missing 
energy on the bond is the excitation energy, whose prob- 
ability distribution, however, is nontrivial. The knowl- 
edge of this distribution would allow to calculate the fi- 
nite size scaling of the slope of the specific heat at Td = 0. 
Due to the simplicity of the boundary droplets the dis- 
tribution could be computed with the dimer algorithm. 
The conditional probability that a bond is occupied given 
its random energy is just the probability p(ij) intro- 
duced above in the calculation of the internal energy, see 



2tt 2 [2tt 2 - (W6)(tt 2 + 4GQ] T$ - 8 [2tt 2 - (&/*>) fr 2 + 16GQ] (£ d /b)T d 

wn+m/bf 2 
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Sec. lIVBl It can be easily obtained from the dimer algo- 
rithm. The distribution of the boundary droplet energies 
is then given by p(ij)p(occ) /p(E) withp(occ), p{E) being 
the probabilities for the occupation and energy = E 
of a bond, respectively. 




linear system size L 

FIG. 14: Scaling of the slope of the dimer specific heat Cd at 
Td = with the system size. 

The smallest droplet excitation in the bulk is likewise 
easily identified as the rotation of a plaquette that con- 
sists of two opposite dimers. The probability distribu- 
tion for the energy difference of the two configurations 
is, however, not easily obtained. It is complicated by 
the condition that the two dimers before the flip must 
be part of the groundstate configuration. In the simula- 
tions, it should be stressed, statistics of droplet energies 
can in principle be measured systematically in the follow- 
ing straightforward procedure. For a given disorder con- 
figuration the groundstate dimer covering is determined. 
Then the energy of one arbitrary occupied bond in the 
bulk is set to infinity and the new groundstate is deter- 
mined. It will not contain the bond with infinite energy 
and hence have higher energy than the original ground- 
state. The energy difference E together with the diam- 
eter I of the non-overlapping region of the two ground 
states is measured. The statistics of these pairs of val- 
ues for many disorder realizations yield the distribution 
pi(E). Quantitative support of the scaling prediction of 



droplet theory is in reach considering the orders of mag- 
nitude over which the dimer model can be simulated. 



V. SUMMARY AND DISCUSSION 

In this article we have compared recent exact replica 
Bethe ansatz results for the planar line lattice with nu- 
merical simulations of the classical random bond dimer 
model. We found excellent agreement for a large set 
of disorder averaged thermodynamic quantities, namely 
the effective disorder renormalized elastic stiffness, the 
free energy, the internal energy, entropy and the specific 
heat of the dimer model. Characteristics of the disor- 
der distribution of thermodynamic quantities can also be 
obtained both from the Bethe ansatz and from simula- 
tions and their agreement has explicitly been shown for 
the variance and the skewness (third cumulant) of the 
free energy. The comparison thus confirms the Replica 
Bethe Ansatz calculation of Ref.lTlland makes the (1 + 1)- 
dimensional line lattice one of the few glassy systems for 
which the validity of a replica approach (without replica 
symmetry breaking) can be critically tested in detail. In 
the comparison, only one free parameter has been used 
which is the ratio of the short scale regularization lengths 
of the continuous line model and the lattice constant of 
the discrete dimer model. The ratio is found to be consis- 
tently of the order of unity for all studied thermodynamic 
quantities. One term in the single line free energy, first 
given in Ref. 0> could not be confirmed by the simula- 
tions, in consistency with numerical data of Ref. [3(]. It 
would be interesting to relate the droplet excitations of 
the dimer or line system to the excitations of the related 
SU(n) fermi gas in the limit n — > 0. In view of the man- 
ifold links of dimer covering models to condensed matter 
systems, especially of spin systematise the link to the 
exactly solved disordered line lattice model might prove 
useful in future applications. 



Acknowledgments 

Interesting discussions with M. Kardar and S. Scheidl 
are gratefully acknowledged. This research was sup- 
ported by the Deutsche Forschungsgemeinschaft (DFG) 
through an Emmy Noether grant (S.B. and T.E.) and by 
NSF Grant DMR-0049176 (A.T. and C.Z.). 



1 K. Binder and A. Young, Rev. Mod. Phys. 58, 801 (1986). 

2 T. Nattermann and S. Scheidl, Adv. Phys. 49, 607 (2000). 

3 Mesoscopic Quantum Physics, edited by E. Akkermans et 
al. (North- Holland, Amsterdam, 1994). 

4 D. S. Fisher, Phys. Rev. Lett. 56, 1964 (1986). 

5 M. Mezard and G. Parisi, J. Physique I 1, 809 (1991). 

6 D. S. Fisher and D. A. Huse, Phys. Rev. B 38, 386 (1988). 



7 C. Zeng, A. A. Middleton, and Y. Shapir, Phys. Rev. Lett. 
77, 3204 (1996). 

8 C. Zeng, P. L. Leath, and T. Hwa, Phys. Rev. Lett. 83, 
4860 (1999). 

9 C. A. Bolle, V. Aksyuk, F. Pardo, P. L. Gammel, E. Zeldov, 
E. Bucher, R. Boie, D. J. Bishop, and D. R. Nelson, Nature 
(London) 399, 43 (1999). 



16 



10 M. Kardar, Nucl. Phys. B 290, 582 (1987). 

11 T. Emig and M. Kardar, Nucl. Phys. B 604, 479 (2001). 

12 T. Emig and M. Kardar, Phys. Rev. Lett. 85, 2176 (2000). 

13 T. Emig and S. Bogner, Phys. Rev. Lett. 90, 185701 
(2003). 

14 N. Elkies, G. Kuperberg, M. Larsen, and J. Propp, J. Al- 
gebr. Comb. 1, 111 (1992). 

15 N. Elkies, G. Kuperberg, M. Larsen, and J. Propp, J. Al- 
gebr. Comb. 1, 219 (1992). 

16 P. W. Kasteleyn, Physica (Utrecht) 27, 1209 (1961). 

17 C. L. Henley, J. Stat. Phys. 89, 483 (1997). 

18 L. D. Landau and E. M. Lifshitz, Elasticity Theory (Perg- 
amon Press, Oxford, 1969). 

19 D. R. Nelson and H. S. Seung, Phys. Rev. B 39, 9153 
(1989). 

20 V. L. Pokrovsky and A. L. Talapov, Phys. Rev. Lett. 42, 
65 (1979). 

21 P. G. de Gennes, J. Chem. Phys. 48, 2257 (1968). 

22 B. Sutherland, Phys. Rev. Lett. 20, 98 (1968). 

23 M. Takahashi, Prog. Theor. Phys. 44, 899 (1970). 

24 M. Kardar and D. R. Nelson, Phys. Rev. Lett. 55, 1157 



(1985). 

25 J.-P. Bouchaud and H. Orland, J. Stat. Phys. 61, 877 
(1990). 

26 E. Brunet and B. Derrida, Phys. Rev. E 61, 200 (1987). 

27 E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963). 

28 N. N. Bogoliubov, J. Phys. USSR 11, 23 (1947). 

29 S. E. Korshunov and V. E. Dotsenko, J. Phys. A 31, 2591 
(1998). 

30 J. Krug, P. Meakin, and T. Halpin-Healey, Phys. Rev. A 
45, 638 (1992). 

31 W. Press, B. Flanery, S. Teukolsky, and W. Velterling, 
Numerical Recipes in C, 2nd ed. (Cambridge Univ. Press, 
1993). 

32 D. S. Fisher and D. A. Huse, J. Phys. A 20L, 1005 (1987). 

33 P. W. Anderson, B. I. Halperin, and C. M. Varma, Phil. 
Mag. 8, 1 (1972). 

34 M. E. Fisher, J. Math. Phys. 7, 1776 (1966). 

35 M. Kac and J. C. Ward, Phys. Rev. 88, 1332 (1952). 

36 R. Moessner and S. L. Sondhi, Prog. Theor. Phys. Supp. 
145, 37 (2002). 



